We work in the context of projector quantum Monte Carlo (PQMC) simulations,^[13,14] where, given an essentially arbitrary initial state |Ψ₀⟩ (often called a “trial state”, though the term is somewhat misleading), the ground state can be found by time evolution in imaginary time,

Although in some cases one would like to reach the true ground state, in studies of continuous quantum critical points it is well known that this is not necessary. Instead, if the dynamic exponent z is known, one can choose β = aL^z, where z is the dynamic critical exponent, with an arbitrary proportionality factor a. This scaling of β removes the dependence on β from the scaling function, and one can then study finite-size scaling of computed quantities only in the spatial size L. This is often a more practical (less demanding of computer resources) alternative to eliminating the β dependence by effectively taking the limit β → ∞ when studying the true ground state.

In analogy with T > 0 simulations of criticality, we will here carry out PQMC calculations with the projection time scaled as τ = aL^z for a quantum-critical system. This has certainly been done before, under the intuitively clear notion that the temporal boundary conditions should not matter and one can proceed exactly as in T > 0 calculations; see, e.g., the recent work in Ref. [18]. However, the freedom to choose the trial state is an aspect of the problem not present in T > 0 simulations, where one always has periodic imaginary time boundaries. To our knowledge, systematic studies of the role of the trial state and how the typicality emerges irrespective of it and the scale factor in the projection time has not been studied systematically before. We will here demonstrate this independence of the asymptotic scaling behaviors on a as well as the trial state |Ψ(0)⟩, by choosing a range of different trial states and factors. Our results in all cases demonstrate quantum-critical typicality of imaginary-time evolved states. On a practical level, this can save time in PQMC simulations, as the simulation time typically scales linearly with τ and reaching the true ground state to within small error bars may require very large τ and detailed tests to check for convergence. Instead scaling τ as aL with a small factor a can then save significant computer time. However, as one may expect, if a is too small, very large system sizes are required for the asymptotic critical behavior to set in. We will observe how the cross-over depends on a and the trial state and also discuss possible elimination of the leading scaling corrections by identifying an optimal value of a. We will also make other interesting observations on the role of the trial state.

To validate the above typicality hypothesis, we use the S = 1/2 Heisenberg model on a symmetric bilayer,^[19–22] illustrated in Fig. 1, as a concrete example. The Hamiltonian of the model is given by

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) The bilayer Heisenberg model with two different exchange constants between nearest-neighbor S = 1/2 spins; J1 within the individual layers and J2 between the layers.

The bilayer model realizes a quantum phase transition from Néel order for small g = J₂/J₁ (e.g., for g = 0 the system consists of two decoupled two-dimensional (2D) Heisenberg layers, for which the long-range Néel order is well understood and quantified^[23]) to a quantum paramagnet for large values. For g → ∞ the ground state is simply a product of singlets on the bond between the layers and there is a spin gap Δ = J₂. This gap closes at g_c, and long-range Néel order forms continuously for g < g_c. By symmetry, this T = 0 transition in 2+1 dimensions (two space dimensions and one time dimension) should belong to the universality class of the finite-T transition of the three-dimensional (3D) classical Heisenberg model, or O(3) model,^[24–26] provided that the Lorentz invariance is emergent when L → ∞ (i.e., the dynamic exponent takes the value z = 1).

The Néel–paramagnetic transition has been studied intensely, with the bilayer model and several other cases of dimerization; for a review, see Ref. [16]. Among the more precise studies of the bilayer, in Ref. [27] SSE calculations were used to determine the critical coupling ratio as g_c = 2.5220(1) (where here and henceforth the number within parenthesis indicates the statistical error of the preceding digit) and the correlation length exponent was found to be ν = 0.7106(9), which is in agreement with the 3D classical Heisenberg exponent ν = 0.7112(5) obtained in a high-precision study of the classical 3D ϕ⁴ model.^[28] This exponent is also consistent with that of the 2D columnar dimerized quantum Heisenberg model.^[29] Initial disagreement with O(3) universality for the case of staggered dimers^[30] have now been attributed to strong corrections to scaling.^[31–33] Thus, there is little doubt that all these dimerized models belong to the same standard O(3) universality class. Here our goal is not to reconfirm this or to obtain more precise exponents, but to convincingly demonstrate that typicality at the quantum critical fluctuations holds, in the precise sense that the same O(3) exponents are produced asymptotically (for sufficiently large system size) with different prefactors a ≤ 1 in the scaling of τ with L in PQMC simulations and with a set of completely different trial states in Eq. (3). We only assume that the dynamical exponent z = 1, and extract the other exponents using finite-size scaling of several physical observables.

In Section 2 we review the PQMC method for quantum spin systems in the valence bond (VB) basis, to make clear the role of the boundary conditions of the time evolution. We also introduce the four trial states used in the simulations and define the physical observables and their corresponding PQMC estimators that we use to study the critical fluctuations. In particular, we discuss the estimator of the spin stiffness, which so far has only been derived within T > 0 QMC methods but for which we here present a simple generalization for PQMC calculations. In Section 3 we discuss the finite-size scaling ansatz within which we analyze our data. We then determine the critical coupling of the bilayer model and extract its universal exponents at criticality. We also discuss atypical and non-asymptotic properties that originate from the initial trial states and finite projecting time. We summarize the results and further discuss them in Section 4.

An antiferromagnetic Heisenberg Hamiltonian can be written as a sum −∑_ijJ_ijP_ij of singlet projectors,

By sampling the spin configurations corresponding to a given singlet tiling of the trial state, one can formulate a PQMC algorithm that is very similar to the T > 0 SSE method running at β = 2τ, including very efficient loop updates.^[35] Let us compare the two approaches. In the SSE method, e^−βH is Taylor expanded up to all contributing orders n, with the maximum contributions of order n = −β⟨H⟩ ∝ βN. The expansion order n is sampled in the simulation according to the total weight of contributions from that order. Often, for practical reasons, a self-selected upper bound n_max is imposed, so that the sampling scheme can be formulated with a fixed number of operators in the operator strings corresponding to the evolution of traced-over states |α⟩ in the chosen basis (normally the basis of spins). Then n of these n_max operators are selected among the terms H_b of the Hamiltonian and the remaining n_max − n entries are “fill-in” identity operators.

Essentially, going from the SSE to the PQMC method corresponds to opening up the periodic time boundaries arising from the trace operation, formally replacing the trace ∑_α ⟨α| ··· |α⟩ by a sum ∑_αL,αR c_Lc_R⟩ α_L|···|α_R⟩ corresponding the projection out of the trial state. In the valence bond basis, this results in “sealing” open loop segments at the time boundaries with VBs. Then, as in the SSE case, a space-time spin configuration can be fully decomposed into loops that can be flipped independently of each other. One can show that all signs cancel out in the overall weights for the PQMC configurations, because of the bipartiteness of the system.

The spin and valence bond representations of a configuration are illustrated by an example in Fig. 2, including the concept of loop updates in the spin representation, Fig. 2(a). The VBs in the trial state can be updated by simple reconfigurations of pairs of bonds, to maintain the A–B connectivity with anti-parallel spins on each bond, using the weights w_r in standard Metropolis acceptance probability. One can also formulate loop updates of the trial state,^[35] but this is more useful in variational calculations than in PQMC.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. (color online) Two different representations used for the configuration space in PQMC simulation in the VB basis, exemplified with a Heisenberg chain of 10 sites. VBs of the sampled trial state are shown at the left and right edges. In panel (a), the linked vertex representation used for loop updates is shown, with solid and open circles indicating up and down spins. Diagonal and off-diagonal operators, and are denoted by open and solid vertical bars, respectively. An example of a loop is shown in red; flipping this loop changes all up spins to down, and vice versa, thereby also changing the operator type (diagonal or off-diagonal) on some of the vertices. The number of operators n and their locations on the lattice are updated separately in “diagonal updates”. In panel (b), the pure VB representation used for collecting operator expectation values is shown for the same configuration. Here there is no distinction between diagonal and off diagonal operators and the bars indicate the full singlet projectors Pij. Propagated left and right VB configurations are shown at the center; they make up the transition graph on which correlation functions are evaluated.

When evaluating operator expectation values, it is normally (depending on the type of operator considered) better to return to the VB only basis, illustrated in Fig. 2(b), which corresponds to summing over all spin configurations that are compatible with the VBs in the trial state and the lattice locations of the operators in the string — here also all the operators turn into the full singlet projectors P_ij, instead of their individual diagonal and off-diagonal terms in the spin basis.

In principle, we do not have to sample VBs of the trial state; we can just use a fixed configuration of the VBs and that will still have overlap with the ground state and converge to the same. However, one can also construct good, practically workable variational states and this can improve the convergence properties.^[35] Even without optimizing, one can write down simple translationally invariant states to further restrict the simulations to total momentum k = 0. In T > 0 simulations, states with all k and all total spin values S are included, and to reach the ground state T has to be well below the smallest gap from the k = 0, S = 0 ground state. PQMC simulations are restricted by construction to S = 0 and k = 0, and, thus, most of the low-energy states are excluded from the outset and do not have to be projected away. Here our purpose is not to reach the ground state perfectly, and we will test the typicality hypothesis both with k = 0 trial states and with simple “frozen” VB states that do not conserve k.

The amplitude-product states proposed by Liang et al.^[13,34] are good k = 0 variational states to describe Néel ordered, critical, and quantum paramagnetic systems. The wave-function coefficients are of the form

As the third trial state we choose a single VB configuration,

The last trial state is again a simple VB configuration, but, unlike |Ψ₃⟩, it breaks the translational symmetry of the system. We choose a columnar arrangement of the VBs in each layer;

The energy expectation values E = ⟨H⟩/N of the above trial states can be evaluated by sampling the VBs of |Ψ₁⟩ and |Ψ₂⟩, and for |Ψ₃⟩ and |Ψ₄⟩ exact calculations are trivial. For L = 16, as an example, the results at the estimated critical point g = 2.5222 (see further below) are E₁ = −1.841(1), E₂ = −1.829(1), E₃ = −1.7611, and E₄ = −1.190275. The unbiased, sufficiently projected energy is E = −1.94225(2). Obviously, the last two trial states are far away from good variational states describing the criticality of the current model. Even |Ψ₁⟩ and |Ψ₂⟩ are not very good variational states, and much better ones can in principle be obtained by optimizing the bond amplitudes. However, our purpose here is not to optimize the states and the simulations, but to demonstrate that typical critical fluctuations emerge out of arbitrary trial states with projection time τ ∝ L. For this purpose the above range of trial states will suffice.

The staggered magnetization, the order parameter, is defined as

The Binder ratio^[40,41] is defined as

The spin stiffness ρ_s characterizes the tendency of ordered spins to adapt in response to a twist imposed on the spins in an ordered state in a direction perpendicular to the ordering vector. In the common QMC simulations at finite temperature, e.g., with the SSE method^[16] or path integrals,^[42] there is a very convenient estimator based on fluctuations of the winding number characterizing the topology of the spin world lines propagated around the space-time periodic system. The stiffness along the α lattice direction is given by

The winding number formula (13) is a consequence of the periodic boundaries in both the spatial and time directions at T > 0. The “cutting open” of the time boundaries in the PQMC method makes this estimator fail at first sight. However, since the current fluctuations within some local time segment Δτ should be independent of β for large β, and these fluctuations are what gives rise to the winding numbers, it is clear that the conservation of the winding number is not very consequential in Eqs. (12) and (13), but is just a byproduct of the periodic time boundaries in combination with the conserved magnetization of the system. Thus, it should be possible to generalize the formula by considering generalized, non-integer winding numbers over a sufficiently large time interval Δτ ≤ 2τ, where we recall that the total length of the system in imaginary time is 2τ. Since the time dimension is not uniform, due to the open boundaries, one can also presume that convergence to the correct value when 2τ, Δτ → ∞ will be faster if Δτ < 2τ and the interval is time-centered.

As a reasonable choice satisfying the above requirements, we take the centered interval with Δτ = τ and calculate the spin stiffness of a d-dimensional system according to

To validate the above formula, we simulate the square lattice Heisenberg model with system size L = 16 using both the SSE and PQMC methods. Results obtained by the two methods are shown versus β and 2τ in Fig. 3. As the ground state is approached for large β and τ, we can see convergence to the same value. We find ρ_s = 0.13239(2) from winding number fluctuations in the SSE simulations and 0.1324(1) from the generalized spin current definition in Eq. (14) in the PQMC simulations. The two estimates agree perfectly. In the case of SSE, we expect asymptotic exponentially fast convergence (though the above result was obtained by extrapolation using a polynomial, which should be fine at the level of the statistical error bars), reflecting the finite-size gap in the spectrum, while in the case of the PQMC calculations with the estimator (14) the convergence appears to be linear in 1/τ. We do not currently have an understanding of this behavior, though clearly it must be related to the fact that the Δτ region over which the current fluctuations are summed have open boundaries and one may expect a correction proportional to the inverse of the length of the boundary. Thus, SSE may still be the preferable way to compute ρ_s, though certainly this test (and others) shows that this important physical quantity can also be reliably obtained in PQMC calculations.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. (color online) Spin stiffness of the 2D Heisenberg model on the square lattice of linear size L = 16 as a function of inverse temperature β in SSE simulations or the projection time 2τ = β in PQMC simulations. The results were obtained using the standard winding number fluctuations according to Eq. (12) and the corresponding generalized spin current formula (14), in the SSE and PQMC simulations, respectively.

For a singular quantity A ∝ δ^κ in the thermodynamic limit, according to the standard finite-size scaling theory,^[43] scaling of the following form is expected close to the critical point g_c:

To locate the critical point, we can treat the scaled A(g,L)L^κ/ν as a dimensionless quantity . By Talyor expanding the scaling function in Eq. (15) and keeping the correction, we have

Alternatively the correlation exponent ν can also be estimated by staying at the size-extrapolated critical point, if this point has been located to sufficient precision. We will use this approach here. First, we calculate the derivative s(g_c, L) of to g at the estimated critical point g_c. This is done by fitting a polynomial f(g) = ag² + bg + c, with three unknown constants a, b, c, to six values of near g_c. Error bars can be estimated by Gaussian noise propagation. Then, according to the following scaling formula,

The Binder ratio Q is dimensionless, which means κ = 0, while for the spin stiffness ρ_s, κ = (d + z − 2) ν, or κ = ν in the present case where d = 2, z = 1.^[46] Therefore Q and ρ_sL are useful observables for locating the critical point and estimating the correlation exponent ν.

For the order parameter m, κ is the exponent β, which leads to the scaling behavior of the squared staggered magnetization,

We next study the typical behavior of the critical fluctuations, by performing PQMC simulations of the bilayer Heisenberg model with the four different trial states defined above and with different prefactors a in τ = aL and using L up to 128. We typically used 10⁵ MC steps to equilibrate the system and 10⁶ for collecting data for the physical quantities of interest. To project out the ground state fully, τ needs to satisfy τ ≫ 1/Δ, with Δ the gap between the ground state and the first excited state “seen” in the calculations, which in the VB basis is the second singlet state. If we can find good scaling properties even significantly away from this limit, it means that a band of low-lying singlets also share the same critical fluctuations as the ground state.

In Fig. 4, we first show results for the squared order parameter evaluated at the critical point, estimated below to be g_c ≈ 2.5222, as a function of the projection time for three system sizes and two trial states. On the scale used in the figure, τ = L gives results almost indistinguishable from the ground state — a close examination (inset of the figure) reveals that there are still some statistically significant differences. With the worst of the trial states |Ψ₄⟩, the results have visibly not converged for τ = L/2, and at τ = L/4 and smaller both trial states give results clearly different from the ground state. Note that one |Ψ₁⟩ has strong Néel that is decays away with increasing τ, while |Ψ₄⟩ has no long-range correlations at all; thus the critical correlations are gradually emergent with increasing τ. Thus, we have a range of different trial states and it is interesting to see if the critical correlations can emerge universally even for small factors a in τ = aL, where the behaviors in Fig. 4 look completely different for the two trial states.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) The squared sublattice magnetization obtained in critical bilayer PQMC simulations with system sizes L = 32,48, and 64, with the trial states |Ψ1⟩ and |Ψ4⟩ (representing the best and worst state as judged by the variational energy). The inset shows the L = 64 results close to τ = L on a more detailed scale.

To analyze the critical point, we begin by considering the time regime where we have almost reached the ground state, using the best trial state in the variational sense, |Ψ₁⟩, and the projection time set to τ = L, i.e., the factor a = 1. The scaled spin stiffness ρ_sL and the Binder ratio Q are shown versus the coupling ratio g for various system sizes in Fig. 5. The results do not differ appreciably from SSE results obtained at very low temperatures,^[27] indicating that the projection time here brings us almost to the ground state. We find crossing points between results for system sizes L and 2L using polynomials fitted to the data points. These crossing points extracted for the two different quantities for a large number of size pairs are shown in Fig. 6(a). The drifts of the g^* values obtained from both quantities are monotonic in L, and both of them converge to a common critical point g_c rapidly for large L. All L ≥ 16 points are consistent with the expected power-law, Eq. (18).

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) The scaled spin stiffness ρs L (a) and the Binder ratio Q (b) of the bilayer Heisenberg model with several L graphed versus the coupling ratio g. The PQMC calculations were carried out using τ = L (a = 1).

Fig. 6.

Figure Option ViewDownloadNew Window

Fig. 6. (color online) Scaling of crossing points g*(L) of the size-scaled spin stiffness ρsL (triangles) and the Binder ratio Q (circles). The results were obtained from four different trial states, as indicated by the legends, and two different projection time factors a were in the PQMC simulations; a = 1 in panel (a) and a = 0.25 in panel (b). The largest system size pair is (L, 2L) = (64,128). The solid curves are fits to g* from ρs L crossings (red) and Q (blue) obtained from states projected out of the best trial state, |Ψ1⟩.

A nonlinear fit of g^*(L) from the ρ_sL crossings yields g_c = 2.52222(4) and 1/ν + ω = 1.65(3), with a reasonable reduced goodness-of-fit value χ² = 1.2. Here the result for the exponent combination 1/ν + ω is not very close to the expected O(3) value, 1/ν + ω ≈ 2.2, likely reflecting the role of remaining higher-order corrections. Such still not size-converged “effective exponents” are known to not significantly effect the extrapolated critical point value.^[45] A similar fit of g^* of Q gives g_c = 2.52224(6) and 1/ν + ω = 2.3(1), with reduced χ² = 1.6. Both estimates of critical point agree well with earlier estimate g_c = 2.5220(1) obtained by using SSE QMC,^[27] in which the ground-state properties were obtained by the β doubling approach, but the statistical error is significantly reduced.

We next consider results obtained with the other trial states: |Ψ₂⟩, |Ψ₃⟩, and |Ψ₄⟩. The projection lengths are first all set as τ = L. The (L, 2L) crossing points g^* from ρ_sL and Q are both shown in Fig. 6(a) together with the previous results based on |Ψ₁⟩. We see only small differences between the results from the different trial states and, not surprisingly, all these crossings converge to the common critical coupling g_c. Though it is not apparent from the figure, somewhat larger system sizes are needed to fit the |Ψ_{t > 1}⟩ results to power-law forms than what is the case with |Ψ₁⟩. The latter state also is the best state in the sense of the variational energy. The results of all the fits are listed in Table 1, including the smallest size used in the fit.

Table 1.

Table 1.

Table 1. Results of finite-size analysis with a single power-law correction to the critical point gc for different trial states and different value of the projection factor a. The standard goodness-of-fit per degree of freedom is denoted as χ2. . a = 1 Trial state gc (ρs) χ2 (ρs) Lmin gc (Q) χ2 (Q) Lmin Ψ1 2.52222 (4) 1.2 8 2.52224 (6) 1.6 16 Ψ2 2.5221 (1) 0.5 12 2.5220 (4) 0.6 12 Ψ3 2.5220 (2) 1.0 12 2.5219 (4) 1.2 12 Ψ4 2.5220 (8) 0.5 12 2.5218 (5) 0.6 12 a = 0.5 gc (ρs) χ2 (ρs) gc (Q) χ2 (Q) Ψ1 2.5215 (4) 0.3 12 2.5221 (1) 0.7 12 Ψ2 2.5218 (3) 1.9 12 2.5221 (4) 0.8 12 Ψ3 2.5220 (2) 0.6 12 2.5221 (1) 1.1 12 Ψ4 2.5220 (2) 0.6 12 2.5221 (2) 1.5 12 a = 0.25 gc (ρs) χ2 (ρs) gc (Q) χ2 (Q) Ψ1 2.5213 (2) 0.9 12 2.5213 (4) 1.7 12 Ψ2 2.5222 (2) 0.8 12 2.5222 (1) 1.4 12 Ψ3 2.5220 (8) 1.4 12 2.5223 (3) 0.4 12 Ψ4 2.5224 (5) 2.3 12 2.5212 (3) 1.4 12

Table 1. Results of finite-size analysis with a single power-law correction to the critical point gc for different trial states and different value of the projection factor a. The standard goodness-of-fit per degree of freedom is denoted as χ2. .

Since the above a = 1 simulations deliver results quite close to the ground state for all the trial states, we need to go to smaller a to investigate the emergence of critical typicality in greater detail. Crossing points obtained from the τ = L/4 simulations are shown in Fig. 6(b). Here we can see very significant differences from the τ = 1 data, but all crossing points still flow toward the same critical point. Results of extrapolations are summarized in Table 1.

We next validate that the states projected out from various trial states at τ ∝ L display typical critical fluctuations characterized by the correct critical O(3) exponents. We demonstrate that this universality emerges with increasing system in finite-size scaling for all the different trial states.

First, we show that the correlation length exponent ν can be extracted from the projected states. We extract the exponent according to Eq. (20) from the g derivatives s(g,L) of the curves Q(g) at the estimated critical point, g_c = 2.5222. The derivatives are extracted from a quadratic polynomial fitted to a set of points for g in the neighborhood of g_c. Results are shown in Fig. 7.

Fig. 7.

Figure Option ViewDownloadNew Window

Fig. 7. (color online) The derivatives with respect to g of ρs L (triangles) and Q (circles) at the best estimated critical point gc = 2.5222, graphed versus the system size on logarithmic scales. Three different projection times τ = aL were used; (a) a = 1, (b) a = 0.5, (c) a = 0.125. Different graphing colors refer to results obtained with different trial states. The solid lines with slope corresponding to the known value of the exponent, ν = 0.7112, in Eq. (20) are draw to show the expected large-L behavior.

In the case a = 1, for all four trial states and for data from both ρ_sL and Q, the derivatives reproduce the expected power-law behavior Eq. (20) as shown in Fig. 7(a). Here we do not include the correction term L^−ω, where ω ≈ 0.78 is expected for the universality class, as we find that the prefactor is small and good fits to just the leading power law L^1/ν can be achieved if some of the smaller systems are excluded. Thus, we extract the correlation exponent ν from the data for each trial state by fitting to Eq. (20), starting from system sizes L_min sufficiently large for the quality of the fit to be acceptable. To avoid the systematic error induced by the small deviations of our g_c value from the true critical point, we here only use system sizes up to L = 64. We have estimated that the deviations will then affect the extracted exponents less than the purely statistical errors of the fitting parameters. All estimated ν values for the four trial states are consistent with each other and with the known value of the exponent. The results are listed in Table. 2. We can see good agreement with the correct O(3) exponent in all cases.

Table 2.

Table 2.

Table 2. Exponent ν obtained from the g derivatives of ρsL and Q at the best estimated critical point gc. To minimize the systematic errors originating from the deviation of the g value used from the true critical point, we only used system sizes up to L = 64 in the fits giving the results shown here. The minimum size is indicated in each case. . a = 1 Trial state ν (ρs) χ2 (ρs) Lmin ν (Q) χ2 (Q) Lmin Ψ1 0.705 (7) 1.3 24 0.716 (7) 1.5 16 Ψ2 0.706 (5) 1.4 24 0.718 (7) 1.3 24 Ψ3 0.707 (7) 1.5 24 0.710 (6) 1.1 16 Ψ4 0.713 (7) 1.3 28 0.711 (5) 0.9 16 a = 0.5 ν (ρs) χ2 (ρs) Lmin ν (Q) χ2 (Q) Lmin Ψ1 0.704 (20) 1.5 32 0.709 (20) 0.9 28 Ψ2 0.709 (8) 0.9 20 0.709 (8) 1.4 12 Ψ3 0.707 (9) 0.8 28 0.714 (9) 0.8 24 Ψ4 0.709 (10) 1.3 24 0.708 (9) 1.1 16

Table 2. Exponent ν obtained from the g derivatives of ρsL and Q at the best estimated critical point gc. To minimize the systematic errors originating from the deviation of the g value used from the true critical point, we only used system sizes up to L = 64 in the fits giving the results shown here. The minimum size is indicated in each case. .

For a = 0.5, the data, shown in Fig. 7(b), show more dependence on the trial state, but in all cases the slope takes the expected value for sufficiently large system sizes. For a = 0.125, we can see very significant dependence on the trial state, but here as well the slope eventually crosses over to the correct critical form for large L. The results for a = 0.5 are also listed in Table 2, but for a = 0.125 we did not carry out the analysis in detail because of the small number of points falling in the asymptotic scaling regime. Nevertheless, these tests make clear that there is a cross-over size, which increases with decreasing a and depends on the trial state, above which the critical O(3) scaling is obtained. The non-universal prefactor of the scaling function depends strongly on a and the trial state.

Next, we investigate the exponent η (the anomalous dimension) of the critical correlation function. We use the squared staggered magnetization ⟨m²⟩ at g_c = 2.5222. In this case we include all our large system sizes in the fits, based on an estimation of the effects of the precision of the g_c values. For a = 1, as shown in Fig. 8(a), the square order parameter ⟨m²⟩ scales well according to the expected critical form Eq. (21) with increasing size L, with only a very weak dependence on the trial state. For the case a = 0.25, significant differences in the values of ⟨m²⟩ can be observed for the different trial states, as illustrated in Fig. 8(b).

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) The squared order parameter at the best estimated critical point gc = 2.5222 versus the system size. for (a) a = 1 and (b) a = 0.5. The solid lines are fits to the large-L data.

Nevertheless, each group of data points from the same trial state forms a straight line on the double-log graph. For large enough L the correction to scaling vanishes and within statistical errors the slopes of the four lines are identical and fully consistent with the known value of η. Fitting to the expected finite-size form, again leaving out system sizes smaller than L_min chosen such that the fits are acceptable, we obtain values listed in Table 3 for a = 1 and 0.5. For a = 0.25, the results are also consistent with the O(3) exponent, but the statistical errors are much larger, due to the large L_min needed in this case, and we do not list the results.

Table 3.

Table 3.

Table 3. Exponent η obtained from finite-size scaling of the staggered magnetization ⟨m2⟩ at the estimated critical point gc, using the different trial states and projection times. The O(3) value of the exponent is η = 0.0375(5).[28] . a = 1 Trial state Ψ1 Ψ2 Ψ3 Ψ4 η 0.036(3) 0.032(3) 0.038(6) 0.039(3) χ2 1.6 1.6 1.8 1.3 Lmin 48 24 18 40 a = 0.5 Trial state Ψ1 Ψ2 Ψ3 Ψ4 η 0.026(6) 0.036(1) 0.037(2) 0.036(4) χ2 1.9 1.6 1.9 0.81 Lmin 16 20 20 56

Table 3. Exponent η obtained from finite-size scaling of the staggered magnetization ⟨m2⟩ at the estimated critical point gc, using the different trial states and projection times. The O(3) value of the exponent is η = 0.0375(5).[28] .

We have shown that typical critical fluctuations emerge out of arbitrary trial states when projecting in imaginary time τ ∝ L^z, instead of projecting fully to the ground state of each finite system. However, as mentioned, from the point of view of the path integral picture of the PQMC simulations, we can understand a as time-space aspect ratio, with the initial trial states corresponding to different kinds of boundary conditions. One can also think of this as a sudden quench, where the Hamiltonian is changed at τ = 0 from the one (some times unknown one) for which the trial state is the ground state to the critical Hamiltonian, followed by time evolution with the latter.

While the critical exponents are independent of the system geometry, the critical value of the Binder ratio Q is universal only in the sense that the value is determined by the dimensionality and symmetry of the system, irrespective of the details of interaction and lattice structure, but under the condition that the boundary conditions and the geometry, e.g., aspect ratio is fixed.^[47,48] This implies that changing the temporal boundary conditions and the aspect ratio should affect the critical value of Q. Therefore it is expected that the critical value of Q in the current PQMC simulations will change with different trial states and factor a in the time scaling τ = aL. We investigate this dependence next.

Figure 9(a) shows the Binder ratio at g_c = 2.5222 versus 1/L in simulations with the four different trial states and two aspect ratios. It is clearly seen that Q_c changes dramatically when a changes from 1 to 0.125 for each of the trial states. The value of Q_c can be found by fitting Eq. (19) to the data. However, for a = 1, the results for Q_c obtained with all trial states seem to converge to the same value when L → ∞; for |Ψ₁⟩, Q_c = 1.290(3), with reduced χ² = 0.9 and starting size L_min = 32; for |Ψ₂⟩, Q_c = 1.289(2) (χ² = 1.1, L_min = 20), for |Ψ₃⟩, Q_c = 1.291(3) (χ² = 1.1, L_min = 20); for |Ψ⟩₄, Q_c = 1.289(2) (χ² = 0.5, L_min = 32). In these fits, the exponent ω in all cases is consistent with the known value ω ≈ 0.78, though with relative statistical errors of about 20% typically. The extrapolated values of Q agree well with Q_c ≈ 1.293(3) previously obtained for an “incomplete bilayer” Heisenberg model (where the intra-layer couplings are missing on one layer), but differs slightly from the result for the complete bilayer Heisenberg model, Q_c = 1.2858(3);^[27] most likely this disagreement is due to an underestimated error bar in the previous calculation.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Critical Binder ratios obtained in simulations with various trial states after projection time τ = aL with (a) a = 1 and (b) a = 0.125. The four curves are power-law fits used to extract the infinite-size values. In panel (a) we only show the fit to data from trial state |Ψ1⟩, for the sake of clarity.

In the case of small a = 0.125, it is clearly seen in Fig. 9(b) that the differences in Q_c obtained with different trial state can be drastic. Nevertheless, Q_c = 1.519(3) from |Ψ₃⟩ and Q_c = 1.523(4) from |Ψ₄⟩ are still very similar; perhaps even identical asymptotically (The exponent ω found in the two fits agree well with 0.78). This may seem surprising, since these trial states are quite different and have different variational energies. The trial states are similar in the sense that they are simple product states of singlets on neighboring sites, but in one case the translational symmetry is broken and in one case it is not. One may speculate that the value of the Binder ratio for small a and L → ∞ is related to the entanglement structure of the trial state. This would be very interesting and deserves further study.